ARPA  Order  Number:  1408 


Program  Code  Number:  9F10 


Contract  t umber:  DAHC19-69-C-0032 


Contractor : 

The  City  College  Research  Foundation 

The  City  College  of  the  City  University  of  New 

Convent  Avenue  at  138th  Street 

CD 

New  York,  New  York  10031 

GO 

£J*ftle:  STRESS  WAVE  PROPAGATION  THROUGH  EARTH-WATER  SYSTEMS 

Second  Project  Report 

© 

< 

Report  No.  72-409-1 

February,  1972 

Principal  Investigator:  Professor  Carl  J.  Costantino 

212-621-2228 


D  D  C 

n\Ef3[?JME{ 

MAR  8  1972 


Contract  Dates:  1  September  1969  to  31  August  1972 

Amount  of  Contract:  $76,123. 


lEtsEtrtris 

B 


This  Research  was  supported  by  the  Advanced  Research  Projects  Agency 
of  the  Department  of  Defense  and  was  monitored  by  the  U.S.  Army 
Research  Office  under  Contract  No.  DAHC19-69-C-0032 . 


Reproduced  by 

NATIONAL  TECHNICAL 
INFORMATION  SERVICE 

Springflvltf,  V*.  23151 


DISTRIBUTION  &ATEMKNT  A 

Approved  for  public  ielecoe* 
Distribution  Unlimited 


Table  of  Contents 


Section  Title 


Pare 


1.0 


2.0 


2.1 


3.0 

3.1 

3.2 

3.3 


3.4 


4.0 


5.0 


Introduction  1 
Governing  System  Eauations  5 
Solution  Procedure  ? 
Numerical  Results  10 
One-Dimensional  Elastic  Consolidation  10 
Triaxial  Elastic  Soil  Consolidation  u 
Triaxial  Coulomb-Mohr  Model  2? 
McCormick  Ranch  Sand  Model  3 6 
Summary  50 
References  51 
Appendix  A  52 


1 


1.0  INTRODUCTION 

This  report  is  the  first  semi-annual  report  on  Con¬ 
tract  No.  DAHC  19-69-C-0032  with  the  Advanced  Research  Projects 

Agency  entitled  "Stress-Wave  Propagation  Through  Earth-Water  Sys- 

* 

terns".  The  fundamental  objective  of  this  study  is  to  develop  nu¬ 
merical  techniques  to  treat  the  general  two-dimensional  stress  wave 
propagation  problem  through  nonlinear  earth  materials  including 
the  effects  of  water  flow  through  the  earth  materials. 

Prior  to  the  beginning  of  this  study,  a  numerical  tech¬ 
nique  was  developed  to  treat  the  dynamic  wave  problem  through  arbi¬ 
trary  nonlinear  media '{Ref.  8  and  7)  without  including  the  ef¬ 
fects  of  water  on  the  propagation  process.  This  numerical  approach 
is  based  upon  the  finite  element  method  of  analysis  and  led  to  the 
development  of  a  large  computer  program  (termed  the  SLAM  Code  for 
identification,  the  acronym  standing  for  Stress  Waves  in  Layered 
Arbitrary  Media)  to  treat  either  the  general  axisymmetric  or  plane 
(stress  or  strain)  geometric  configuration.  The  finite  element 
approach  has  been  taken  in  this  development  to  allow  the  user  a 
general  flexibility  in  treating  two  dimensional  problems  of  rather 
complex  geometry  (inclusions,  material  layering,  complex  bounda¬ 
ries,  etc.) 
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advantage  can  be  considered  a  disadvantage  for  these  cases. 

After  the  development  of  SLAM  Code,  various  problems  of 
interest  were  investigated  to  determine  the  effects  of  material 
nonlinearities  on  the  wave  propagation  process  (Ref.  5).  In  gen¬ 
eral,  two  types  of  problems  are  of  interest  when  studying  dynamic 
processes  through  earth  media.  In  the  first  type,  the  half  space 
is  subjected  to  high  intensity  pressure  loadings  caused  by  high 
energy  explosions.  The  resulting  ground  shock  effects  are  highly 
transient  and  are  characterized  by  relatively  short  duration  shock 
waves  of  high  strength.  Of  particular  interest  for  this  problem 
is  the  rate  of  deca};  of  the  shock  front  as  it  moves  through  the 
ground.  Clearly,  nonlinear  properties  of  the  material  signifi¬ 
cantly  influence  the  decay  of  the  shock  strength  since  large  non- 
recoverable  volume  changes  can  decrease  the  peak  pressures  of  the 
shock  front.  In  the  second  problem  type,  the  half -space  is  sub¬ 
jected  to  long  duration  low  intensity  vibratory  type  loadings  as¬ 
sociated  with  earthquake  motions .  Again  nonlinear  properties  and 
volume  changes  of  the  earth  material  significantly  alter  the  char¬ 
acteristics  of  the  motion  histories  sustained  at  the  surface  of 
the  ground. 

In  both  of  these  problems,  the  stress  and  motion  histo¬ 
ries  sustained  at  any  point  in  the  ground  are  significantly  influ¬ 
enced  by  the  nonlineat  characteristics  of  the  materials  and  the 
associated  volume  changes  that  occur.  However,  for  real  problems 
of  interest,  the  soil/rock  media  of  ter.  contains  entrapped  pore 
fluid  which,  due  to  its  relatively  high  stiffness,  will  delay  these 
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volume  changes  from  occurring,  the  amount  of  the  delay  being  con¬ 
trolled  by  the  imperviousness  of  the  earth  material  to  water  flow. 
Thus  to  serioucly  treat  the  time  dependent  response  of  earth  media, 
the  effects  of  pore  water  must  be  suitably  taken  into  account. 

In  the  following  developments,  the  effects  of  pore  water 
arp  included  in  the  finite  element  analysis  of  SLAM  Code  with  the 
eventual  goal  of  treating  the  complete  dynamic  process.  In  the 
first  effort  presented  herein,  only  the  quasi-static  problem  is 
considered;  that  is,  inertial  effects  are  neglected.  The  problems 
of  concern  then  are  limited  to  (time-dependent)  two  dimensional 
consolidation  situations  including  material  nonlinearity  effects. 
This  course  has  been  taken  as  a  first  step  on  the  route  to  the 
complete  development.  The  analyses  and  associated  computer  code 
developments  can  then  be  checked  or  compared  with  known  solutions 
already  available.  At  the  end  of  this  comparison  phase,  the  inclu¬ 
sion  of  inertial  effects  can  then  be  completed  and  solutions  ob¬ 
tained  for  the  various  wave  problems  of  interest. 
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2.0  GOVERNING  SYSTEM  EQUATIONS 

In  Ref.  1,  the  derivation  of  the  system  of  equations 
governing  both  the  equilibrium  of  the  nodes  as  well  as  the  pore 
pressure  seepage  condition  (based  on  the  two-dimensional  Darcy's 
Law)  were  presented  in  detail.  These  relations  can  be  written 
symbolically  in  matrix  form  as 

(a)  Equilibrium: 

:  (j*  j.  jW  + 1  wH")  - 

{ u -Ife-JW +  U-J  w  -  (1) 

(b)  Seepage  Flow; 

{  p]  *  [k  ♦  t\„T  W +  LhI 

(2) 

In  equations  1,  the  vectors  (  F,^  and  are  the  horizontal 

(^-direction)  and  vertical  (w-direction)  forces  applied  at  the 
node  points.  The  matrices  kU0L ,  kUur,  kUju>,  k^^are  the  usual  elastic 
stiffness  matrices  and  are  used  to  compute  the  elastic  components 
of  the  forces  developed  at  the  node  points  due  to  relative  dis¬ 
placements  of  the  element  nodes.  The  vectors  and  \ur\  are  the 
horizontal  and  vertical  displacements  of  the  nodes  of  the  mesh, 
while  the  vector  |tt^  represents  the  excess  pore  pressures  developed 
at  the  node  points.  The  matrices  [  and  then  convert  the 

excess  pore  pressures  developed  at  the  nodes  into  equivalent  node 
point  forces. 

The  forces  at  the  node  points,  {  Fu\  and  [fJ\  ,  have  two 
components ,  namely , 
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where  f  are  the  horizontal  and  vertical  components  of 

any  forces  applied  to  the  nodes  (from  concentrated  loads  or  pres¬ 
sures  applied  to  specific  surfaces  in  the  problem)  .  The  forces 

and  ^F*j  are  the  fictitious  correction  forces  that  are  ap¬ 
plied  to  the  nodes  to  account  for  any  nonlinearities  in  material 
stress-strain  behavior  (or  deviations  from  the  elastic  case) .  For 
completeness  of  this  report,  the  formulation  of  the  matrices  are 

presented  in  Appendix  1. 

Equations  1  and  3  then  represent  the  equilibrium  of  total 
stresses  at  a  point  in  the  half-space.  Equation  2  represents  the 
equation  controlling  the  rate  of  seepage  through  the  body  and  is 
obtained  from  Darcy's  Law.  The  Matrix  [h']  is  dependent  upon  the 
coefficients  of  permeability  of  the  earth  material  (as  well  as 
properties  of  the  finite  element  configuration  used)  and  the  vec¬ 
tors  {$\  and  ^  represent  the  horizontal  and  vertical  node  point 

velocities . 

The  vector  represents  the  rate  of  volume  change  of 

the  fluid  associated  with  each  node  point.  For  incompressible 
fluid,  these  components  are  zero  for  interior  node  points  (all 
fluid  that  flows  into  an  element  must  flow  out) ,  while  for  some 
boundary  node  points  (for  which  the  excess  pore  pressure  is  zero) , 
these  components  indicate  the  volume  of  water  flowing  out  of  the 
nodes.  For  compressible  pore  water  (due  to  the  fluid  compressibility 


itself  or  entrapped  air)  ,  these  components  indicate  the  amount  of 
volume  change  undergone  by  the  fluid  (see  Ref.  1). 

2.1  Solution  Procedure 

At  each  node  point  (except  at  boundary  nodes  where  either 
displacements  and/or  excess  pore  pressures  are  specified) ,  three 
unknowns  must  be  determined  at  any  instant  of  time,  namely  the  two 
node  displacements  (a  and  w)  and  the  excess  pore  pressure  OT) . 

The  solution  is  then  marched  out  in  time  in  a  step  by  step  fashion 
The  integration  procedure  used  to  obtain  the  numerical  results  to 
be  presented  is  based  on  a  simple  linear  velocity  approximation 
during  a  "small"  time  step  or 


where  X i  represents  a  displacement  at  time  i,  x^..,  the  displace 
ment  at  the  previous  time,  x;  and  *•_,  represent  the  corresponding 
velocities  and  At  is  the  time  increment  between  l-i  and  1  .  Solving 
equation  4  for  the  current  velocity 

■x-i  *  -  y-i_(  +•  i  (Xf  (5) 

Substituting  equation  j  into  equation  2,  the  seepage 
equilibrium  equation  can  be  written  as 


where  the  vector  is  defined  as 


a 


(7) 


The  subscript  i  in  equations  6  and  7  represents  the  current  time 
and  (-1  represents  the  previous  time.  At  the  current  time,  then, 
the  system  of  equations  to  be  solved  can  be  written  as 


or 


(9) 


The  matrix  [kJ  of  equation  9  is  symmetric,  and  the  usual  solution 
procedures  can  be  used. 

Considering  an  elastic  porous  material,  at  a  particular 
time,  the  effective  force  vector  of  equation  9  is  known.  The  ap¬ 
plied  loads  are  specified  at  the  current  time  (F^  ,  F^.  )  and  the 
force  component  G  can  be  computed  from  equation  7  since  the  solu¬ 
tion  from  the  previous  time  step  is  known  (obviously,  the  solution 
must  start  from  a  time  when  the  initial  conditions  are  specified) . 
The  unknowns  (  of  equation  9)  can  then  be  Obtained  by,  say,  a 
simple  elimination  technique.  The  solution  at  the  following  time 
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step  can  then  be  obtained  using  the  current  solution  as  input,  etc. 
In  this  fashion,  the  solution  is  marched  out  in  time. 

For  nonlinear  material  behavior,  this  process  has  to  be 
modified  since  the  nonlinear  correction  forces  (F^  ,  )  of  equa¬ 

tion  8  are  also  functions  of  the  current  displacements.  To  over¬ 
come  this  situation,  a  modification  of  the  above  procedure  is  nec¬ 
essary.  An  initial  trial  solution  is  first  obtained  by  using  ap¬ 
proximate  values  for  these  correction  forces  (usually  the  forces 
from  the  previous  time  step) .  An  iteration  procedure  is  then  super¬ 
imposed  at  each  time  step  to  check  the  adequacy  of  the  trial  non¬ 
linear  force  correction  terms. 
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3.0  NUMERICAL  RESULTS 

With  the  developed  computer  program,  numerical  results 
were  generated  for  several  soil  configurations  similar  to  the  usual 
soil  tests,  consolidation  and  triaxial  compression.  The  first  set 
of  data  assumed  elastic  soil  behavior,  since  for  these  problems 
analytic  solutions  are  available  or  can  be  easily  developed  for  com 
parison  purposes.  The  first  nonlinear  soil  model  investigated  made 
use  of  a  Coulomb-Mohr  elastic  plastic  model  based  on  the  concepts 
of  the  theory  of  plasticity. 

Although  this  model  is  often  used,  it  is  not  adequate  for 
modeling  stress-strain  behavior  (except  in  a  crude  sense)  and  would 
be  of  questionable  value  when  studying  pore  pressure  dependent  prob¬ 
lems.  A  more  detailed  soil  model  was  then  investigated  which  ade¬ 
quately  predicts  stress  strain  behavior  of  a  particular  sand  sample 
and  was  developed  by  fitting  the  parameters  of  this  model  to  avail¬ 
able  experimental  data. 

3.1  One-Dimensional  Elastic  Consolidation 

The  first  problem  investigated  was,  naturally,  that  of 
the  classical  one-dimensional  consolidation  of  elastic  material. 

The  analytic  solution  available  for  comparison  is  the  standard 
Terzaghi  solution  (Ref.  2).  The  problem  parameters  chosen  for  the 
investigation  are  shown  in  Figure  1.  The  computed  settlement-time 
history  at  the  top  of  the  soil  surface  is  shown  in  Fig.  2  and  com¬ 
parisons  made  with  the  exact  analytic  solution.  The  excess  pore 
pressures  developed  at  the  bottom  of  the  layer  are  shown  in  Fiaure  3 
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while  those  at  a  point  nearer  the  surface  are  shown  in  Figure  4. 
Again  comparisons  are  made  with  the  exact  solution  and  in  all  cases 
they  show  excellent  agreement.  The  pore  pressure  distribution  at 
various  times  through  the  layer  is  shown  in  Figure  5.  Since  the 
pore  pressure  is  assumed  to  vary  linearly  within  a  given  element 
the  distribution  curves  are  piecewise  linear.  If  in  the  actual 
problem  the  pore  pressure  variation  is  sharp,  smaller  element  sizes 
must  be  used  to  suitably  approximate  the  solution. 

3 . 2  Triaxial  Elastic  Soil  Configuration 

The  second  model  considered  was  the  triaxial  soil  con¬ 
figuration  shown  in  Fig.  6a.  The  soil  model  was  considered  to  be 
elastic  and  a  50  psi  vertical  pressure  appLied  at  the  initial  or 
zero  time.  The  finite  element  model  used  is  shown  in  Fig.  6b  and 
consists  of  28  rectangular  elements  to  represent  the  upper  quarter 
of  the  triaxial  sample.  The  elements  are  thus  axisymmetric  or  ring 
elements . 

To  obtain  the  analytic  solution,  it  was  assumed  that 
strain  conditions  in  the  sample  are  uniform.  The  initial  pore  pres¬ 
sure  developed  in  the  sample  (prior  to  drainage  occurring)  is  found 
from  the  following  analysis.  The  volume  change  per  unit  soil  volume 
for  the  elastic  soil  is 

AV«  7”  l<Tr»ff.+ir,]  U0) 


where  the  barred  stresses  represent,  the  intergranular  stresses,  and 


O  ExOct  Solution 


One  Dimensional  Elastic  Consolidation 


-  Element  Z. 


15 


r n 
O 
<D 
vf» 

W 

e 

v- 


c 

o 


•rl 

+) 

+» 

c 

ns 

as 

T3 

g 

•H 

as 

i — 1 

H 

O 

w 

PI 

f! 

•d 

O 

c 

O 

eg 

o 

+» 

•H 

<d 

■P 

(0 

as 

nJ 

M 

i— i 

d 

W 

in 

(0 

i— l 

as 

—4 

‘w 

P 

C 

& 

O 

•H 

as 

(0 

p 

C 

0 

as 

c 

A 

■5 

(0 

Q 

M 

as 

as 

0 

c 

X 

o 

w 

•rt 


>1 


Element 

^presentation 


'V  "  r-- 

E  =  1000  ps>> 

Vs  o. £5  4 

*/«■»='•*  %-sec 


(a)  Conf igurat ion  of  Triaxial  Model 


Free 

Surface 


T 

T 

’ 

1' 

A 

a 

16 

22 

Z 

9 

16 

23 

3 

10 

17 

24 

4 

M 

Id 

25^ 

5 

1 2 

19 

Q> 

20 

27 

7 

14 

21 

28 

^e*v-  bop&i 


^Element  Number 


^  Roller 
Supports 


(S')  Axiej metric  Finite  Element  Model 


Fig.  6  Elastic  Triaxial  Test  Configuration 


E  and  1>  are  the  elastic  modulus  and  Poisson's  ratio,  respectively. 
Since  no  seepage  occurs  during  the  initial  conditions,  the  volume 
change  is  zero  or 


CTe)  (11) 

In  addition, 

=  cfs-- 

(l2) 

where  p  is  the  excess  pore  pressure  and  <3^-  is  the  vertical  applied 
stress.  Combining  equations  11  and  12  leads  to  the  solution 

5,  -  1 V3 

3  '  (13) 

-p  =  °^/3 

The  initial  compression  of  the  soil  sample  is  simply 

Y'  U+H)  (14) 

The  final  stresses  in  the  soil  system  are  obtained  when 


p  is  zero  (no  pore  pressure)  and 


1 J 

while  the  final  compression  of  the  soil  sample  is 


(16) 


where  L  is  half  the  original  sample  height  (height  of  the  finite 
element  model) .  The  settlement  from  the  initial  condition  to  the 
final  condition  is  governed  by  the  one-dimensional  consolidation 
model  (since  one-dimensional  seepage  occurs  through  the  top  surface 

only)  with  the  modification  that  the  definition  of  the  coefficient 
of  consolidation  is 


c 

'  ,17) 

The  solution  to  the  particular  problem  of  Fig.  6  was  ob¬ 
tained  numerically  using  a  time  increment  of  0.1  seconds.  The  pore 
pressure  distribution  along  the  centerline  elements  is  shown  in 
Fig.  7  together  with  comparisons  with  the  analytic  solution.  As 
can  be  seen,  the  comparisons  are  excellent,  except  during  the  early 
part  of  the  solution.  In  an  attempt  to  uncover  the  cause  of  the 
discrepancies,  the  same  problem  was  investigated  with  differing 
time  increments,  and  the  results  are  shown  m  Fig.  8.  As  may  be 
noted  by  comparing  Figs.  7  and  8,  the  early  time  oscillations  found 
for  the  top  element  (Element  1)  are  related  to  the  time  step.  As 
the  time  step  is  decreased,  the  oscillations  disappear.  A  compar¬ 
ison  with  the  exact  solution  shows  that  the  computed  solution  is 
slightly  lower  and  this  can  be  attributed  to  the  fact  that  the 


pore 


Fig.  7  Elastic  Triaxial  Test,  Excess  Pore  Pressure  Near  Centerline 
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pressure  profile  is  assumed  to  be  linear  across  the  element  while 
the  actual  pressure  profile  is  curved,  pacticulirly  at  the  early 
times . 

The  comparison  with  the  middle  element  (Element  4)  is  not 
as  clear  cut,  however.  As  may  be  noted  from  Fig.  8,  this  pore  pres 
sure  shows  an  initial  increase  in  pore  pressure  before  the  antici¬ 
pated  decay  occurs,  and  this  increase  is  independent  of  time  incre¬ 
ment  of  the  integration.  Since  this  phenomenon  did  not  occur  in 
the  elastic  plane  problem  discussed  previously,  it  must  be  concluded 
that  this  variation  is  concerned  with  the  coarseness  of  the  finite 
element  mesh  in  the  radial  direction  for  this  axisymmetric  problem. 
No  further  numerical  studies  have  been  conducted  on  this  problem 
as  yet,  however. 

3. 3  Triaxial  Coulomb-Mohr  Model 

The  first  triaxial  problem  including  nonlinear  material 
properties  that  was  investigated  was  the  same  model  shown  in  Fig.  6 
but  with  nonlinear  properties  described  by  the  Coulomb-hohr  yield 
condition  (Ref.  3).  For  stresses  within  the  yield  surface,  the 
soil  is  assumed  to  behave  elastically,  where  the  yield  surface  is 
defined  by 


dL  J,  ♦  [jj  *  k 


(18) 


For  the  axisymmetric  stress  condition  of  interest  for  this  problem, 


J.  -  <Tr  +  (T0  CTj 


2  j 

(19) 


where  the  bar  again  indicates  intergranular  stresses.  The  coeffi¬ 
cients  (<X,  k)  are  related  to  the  usual  strength  parameters  obtained 
from  a  tnaxial  test  series,  4>  ,  the  angle  of  internal  friction,  and 
c,  the  cohesion,  by 


e*  » 


2_ 

3 


am  4> 
C3-Sir\4>') 


k 


(20) 


For  stresses  on  the  yield  surface,  plastic  strain  components  are 
determined  from  the  usual  normality  principal. 

Prior  to  investigating  this  problem  numerically,  the  ana¬ 
lytic  solution  for  the  initial  stress  condition  was  obtained  (no 
drainage  allowed).  As  the  vertical  stress  is  slowly  increased,  the 
soil  behaves  elastically  and  the  previous  solution  applies.  Sub¬ 
stituting  equations  12  and  13  into  equations  19  yields 


r,  *  o 


For  elastic  yielding  to  begin,  the  critical  vertical  stress  must 


reach  the  value 


qj  *  fe. 
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(21) 


For  applied  stresses  larger  than  this  critical  value,  plastic  flow 
must  be  accounted  for,  making  use  of  the  normality  relation 
(Ref.  3),  which  for  this  problem  becomes 


ay) 


(22) 


•f 


where  (  £r  ,  )  are  the  radial  and  vertical  components  of  the  plas¬ 

tic  strain  rate  vector.  The  plastic  volume  change  is 


(23) 


P 

where  6^  is  the  total  plastic  vertical  strain,  while  the  elastic 
volume  change  is 


4ve.(':r  )(v^r) 


(24) 


Knowing  that  the  total  volume  change  is  zero  (no  drainage  out  of 
the  sample  is  allowed) ,  the  solution  can  be  readily  obtained  for 
any  applied  stresses  greater  than  the  critical,  or 


or 


(25) 


The  results  for  a  particular  undrained  case  are  shown  in 
Fig.  9.  The  vertical  pressure  is  applied  "slowly"  with  a  rise  time 
of  50  seconds  until  it  reaches  a  peak  pressure  of  50  psi.  The  par¬ 
ticular  properties  of  the  soil  chosen  were 


E  -  10 OO  psi 

1)  *  0.25 

e  »  2o.8  •pa*.' 

<t>  *  3 c6 

For  this  condition  the  critical  vertical  stress  is  reached  when  CT„ 
is  43.4  psi  and  the  corresponding  pore  pressure  is  14.45  psi.  As 
the  vertical  stress  is  increased  to  50  psi,  plastic  flow  takes  place 
(along  with  plastic  volume  expansion)  and  the  pore  pressure  reduces 
to  11.1  psi.  Five  computer  runs  were  made  for  this  problem  using 
different  time  steps  as  seen  in  Fig.  9. 


In  each  case,  the  nonlinear 
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Fig.  9  Triaxial  Compression,  Coulomb-Mohr  Material,  Undrained  Test 
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correction  forces  in  the  equilibrium  equations  were  taken  as  the 
value  computed  during  the  previous  time  step.  As  can  be  noted, 
the  smaller  the  time  step,  the  better  the  approximation,  as  expected 
As  an  alternate  to  this  procedure,  the  nonlinear  correction  forces 
in  a  given  time  step  can  be  recomputed  by  iteration  (obtain  a  trial 
solution,  computed  correction  force,  obtain  new  solution,  etc.). 

For  this  problem  of  proportional  loading,  this  procedure  is  equiva¬ 
lent  to  using  smaller  time  steps  without  iteration  during  each  time 
step. 

After  the  final  equilibrium  condition  is  reached  under  no 
drainage  conditions,  the  drained  situation  can  be  achieved  by  let¬ 
ting  the  pore  pressure  decrease  to  zero  by  allowing  drainage  through 
the  top  and  bottom  surfaces  of  the  soil  sample.  It  can  be  shown 
that  for  this  soil  model,  the  decay  of  the  pore  pressure  will  occur 
elastically;  that  is,  the  intergranular  stress  state  will  move  off 
the  yield  surface  as  the  pore  pressure  decreases,  so  that  the  decay 
rate  will  be  as  described  in  the  previous  elastic  triaxial  solution. 

The  solutions  for  these  cases  are  shown  in  Fig.  lo  where 
the  vertical  intergranular  stress  is  plotted  as  a  function  of  the 
total  vertical,  strain  for  various  values  of  the  cohesion  and  for  a 
fixed  value  of  the  friction  angle  of  30*.  If  the  cohesion  is  24.0 
psi  or  greater,  the  soil  sample  always  remains  elastic.  The  ini¬ 
tial  stress  state  when  a  vertical  stress  of  50  psi  is  applied  and 
no  drainage  is  allowed  is  =  33.3  psi  and  p  >  16.7  psi.  When 
drainage  is  then  allowed,  the  pore  pressure  decreases  to  zero,  and 
the  vertical  intergranular  stress  increases  to  the  applied  stress 
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of  50  psi.  The  final  state  is  the  same  as  would  occur  if  the  sample 
had  been  tested  dry  (no  pore  pressure).  If,  however,  the  cohesion 
is  lower,  the  initial  undrained  state  causes  plastic  flow  to  occur 
with  the  drained  condition  occurring  elastically,  as  shown.  For  a 
value  of  cohesion  equal  to  14.4  psi,  the  initial  undrained  state 
occurs  with  no  excess  pore  pressure,  and  the  drained  state  is  the 
same  as  the  undrained  state.  For  values  of  cohesion  less  than  14.4 
psi,  equilibrium  under  the  applied  loads  cannot  be  maintained.  It 
should  be  pointed  out  that  for  values  of  cohesion  between  14.4  and 
24.0  psi,  the  dry  test  will  show  no  plastic  flow,  while  the  undrained- 
drained  sequence  will  yield  plastic  strains. 

it  is  clear  then  that  even  for  this  relatively  simple 
soil  model,  the  stress-strain  behavior  between  saturated  and  unsat- 
uiated  soil  samples  will  be  different  and  will  be  influenced  by 
the  rate  of  loading  (as  compared  with  the  rate  of  pore  pressure 
decay).  To  investigate  this  analytic  solution  further,  the  pre¬ 
vious  solution  was  nondimensionalized  in  the  following  fashion. 


Non  dimensional  parameters  are  defined  as 


fa 


(26) 


The  upper  and  lower  limits  of  cohesion  for  which  a  nonlinear  solu¬ 
tion  (stable  plastic  strains  will  occur)  can  be  obtained  for  the 
undrained  case  are 
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(27) 


For  any  value  of  cohesion  between  these  limits,  the  solution 
yields 
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After  this  initial  solution  occurs,  the  addition  vertical  strain 
that  will  develop  as  the  excess  pore  pressure  is  allowed  to  decay 
to  zero  is 


(29) 


so  that  the  final  strain  is  the  sum  of  the  strains  from  equations 
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28  and  29.  The  solutions  for  several  parameter  variations  are  shown 
in  Pigs.  11  through  14.  In  Fig.  n,  the  nondimensional  vertical 
intergranular  stress  is  plotted  as  a  function  of  the  ratio 

where  is  the  vertical  strain  that  would  occur  in  the  dry  state 
and  is  simply 


As  may  be  noted,  the  difference  in  limiting  values  of  cohesion  for 
this  problem  is  relatively  small,  but  the  influence  on  the  final 
strain  is  large  (ratio  of  6.25).  Curves  are  shown  for  four  equally 
spaced  values  of  cohesion  between  the  limiting  values. 

The  same  solution  is  shown  in  Pig.  12,  except  that  the 
friction  angle  was  increased  from  5°  to  30°.  As  may  be  noted,  the 
final  strains  are  much  lower  than  those  of  Fig.  u,  and  the  associ¬ 
ated  plastic  strains  occurring  during  the  initial  undrained  state 
are  much  smaller.  This  is  due  to  the  fact  that  for  the  higher  fric¬ 
tion  angle  the  plastic  volume  expansion  is  larger  than  for  the 
smaller  friction  angle  causing  the  excess  pore  pressure  to  decay 
more  rapidly  as  plastic  strains  develop.  Fig.  13  shows  the  same 
results  for  a  still  larger  friction  angle  of  45°,  again  showing  a 
smaller  difference  in  final  strains. 

The  results  for  a  different  value  of  Poisson's  ratio 

°-25)  are  shown  in  Pi9-  H  for  a  friction  angle  of  30°.  As 
can  be  noted,  the  behavior  is  essentially  different  than  that  of 
Fig.  12.  This  is  due  to  the  fact  that  the  elastic  volume  change 
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Fig.  11  Triaxial  Compression  Test,  Coulomb-Mohr  Material,  Nondimensional  Solution 
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during  the  initial  loading  decreases  as  Poisson's  ratio  increases. 

3.4  McCormick  Ranch  Sand  Model 

It  is,  of  course,  well  known  that  the  relatively  simpli¬ 
fied  constitutive  models,  such  as  the  Coulomb-Mohr  model,  can  only 
crudely  approximate  the  stress-strain  behavior  of  real  soils.  In 
order  to  properly  take  into  account  the  influence  of  pore  fluid  on 
soil  response,  more  realistic  models  must  be  developed.  An  example 
of  such  a  model  was  presented  in  Ref.  4  wherein  the  parameters  of 
the  model  were  chosen  to  match  (as  closely  as  possible)  available 
experimental  data  on  a  particular  sand  sample,  known  as  McCormick 
Ranch  Sand.  A  rather  extensive  series  of  triaxial,  uniaxial,  and 
hydrostatic  compression  tests  were  conducted  and  an  attempt  was 
made  to  fit  the  analytic  model  so  as  to  reproduce  the  available  data. 

It  was  found  that  for  the  particular  parameters  chosen 
the  stress-strain  curve  during  the  initial  load-unload  cycle  could 
be  adequately  reproduced  for  the  triaxial  compression  test  (over  a 
wide  range  of  confining  pressures)  and  for  the  uniaxial  compression 
test.  The  soil  model,  however,  was  significantly  stiff er  under 
hydrostatic  compression  (although  the  shape  of  the  load-unload  curve 
was  the  same)  than  the  experimental  data. 

The  model  is  based  on  the  following  analysis.  The  hydro¬ 
static  and  deviatoric  stress-strain  components  are  related  by 

Sij  -  1G  €ij 


-p  =  3K  e 


(31) 
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where  ${j  =  deviator ic  stress  tensor 
=  deviatoric  strain  tensor 
“jo  *  hydrostatic  pressure 
g  =  volumetric  strain  *  ^  -f.  ^  +■  ) 

and  are  related  to  the  total  stress-strain  components  by 


•  <tj  -  e  Sij 
Sij  =  %•  -  f  &{j 


(32) 


where  (  Ojj  ,  )  are  the  total  stress-strain  tensors  and  is  the 

Kronecker  delta.  The  dots  in  equations  31  and  32  indicate  the  cor¬ 
responding  rates.  The  parameters  K  and  G  represent  the  bulk  and 
shear  moduli,  respectively,  and  are  taken  as  functions  of  stress 
history. 

The  form  used  for  the  bulk  modulus  is: 


loading:  Kt  =  K0  +  K^e  +  for  k>0  (33) 

unloading:  Ku  =  K6yw  + 

where  the  parameters  K0,  ,  KO0i,  KilA  are  parameters  found  by 

fitting  the  experimental  data.  In  equation  33,  volume  compression 
is  assumed  to  be  positive.  The  corresponding  form  used  for  the 
shear  modulus  is 
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where  -J°c  is  a  critical  hydrostatic  pressure  (positive  in  compres¬ 
sion  and  Jj  is  the  second  invariant  of  the  deviatoric  stresses 
(equation  19) . 

To  match  the  specific  test  results  for  the  sand  sample, 
the  following  parameters  were  found  to  best  reproduce  all  the  data 
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The  stress-strain  behavior  under  uniaxial  compression  as 


shown  in  Fig.  15  under  both  initial  loading  conditions  as  well  as 
strain  load/unload  cycling.  As  may  be  noted,  the  stress-strain 
response  exhibits  the  characteristic  stiffening  effect  as  well  as 
the  nonrecoverable  behavior  under  load  cycling.  The  pressure  ranges 
shown  are  higher  than  normally  used  but  suitable  modification  of 
the  data  input  would  convert  this  typical  response  to  lower  stress 

ranges  of  interest. 

The  behavior  under  triaxial  compression  is  presented  in 
Figs.  16  to  19  and  again  exhibits  much  of  the  characteristics  antic¬ 
ipated  for  a  sand  sample.  During  the  load/unload  cycling,  the 
model  can  be  further  improved  to  reproduce  test  data  by  modifying 
the  shear  nr  lulus  formulation  under  reload  conditions  to  better 

match  strain  behavior  with  constant  load  cycling. 

The  previous  data  were  obtained  for  the  Ranch  Sand  model 

in  the  dry  condition.  To  determine  the  behavior  with  pore  fluid, 
similar  problems  were  investigated  including  load  cycling  effects. 
In  Fig.  20,  the  triaxial  response  is  presented  for  a  consolidated/ 
undrained  experiment  with  load  cycling  in  the  vertical  direction 
corresponding  to  the  load  cycles  shown  in  Fig.  18  for  the  dry  sam¬ 
ple.  In  both  cases,  lateral  or  confining  stresses  were  maintained 
constant.  As  can  be  seen  in  Fig.  20,  the  effect  of  pore  pressure 
is  to  decrease  the  axial  s -rain  increment  between  load  cycles. 

That  is,  in  the  undrained  state,  the  soil  model  "shakes  down"  to 
effectively  a  linear  model,  although  strong  nonlinear  behavior 
again  takes  effect  as  the  apolied  load  is  finally  increased  beyond 
the  load  cyclina  reaime. 
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Similar  behavior  is  shown  in  Fig.  21  where  the  applied 
vertical  load  is  cycled  through  the  complete  load  range  from  0  to 
300  psi.  This  test  corresponds  to  the  dry  triaxial  test  shown  in 
Fig.  19.  Again,  it  may  be  noted  that  within  a  load  cycle,  pore 
pressure  effects  cause  the  stress-strain  behavior  to  "shake  down" 
to  an  effective  elastic  state.  Of  course  this  type  of  response 
can  be  modified  by  changing  the  definition  of  the  reload  shear  mod¬ 
ulus  as  defined  by  equation  34.  A  plot  of  the  invariants  of  effec¬ 
tive  stresses  during  the  loading  cycle  for  the  triaxial  tests  is 
shown  in  Fig.  22,  for  both  the  consolidated  undrained  and  drained 
tests.  As  may  be  noted,  J,  is  constant  during  the  undrained  test 
indicating  that  the  bulk  modulus  (equation  33)  is  constant  with  this 
model.  Therefore  the  cycling  resr-,se  will  be  completely  dependent 
upon  the  variation  in  the  deviatoric  response,  or  the  shear  modulus 
behavior.  The  cycling  response  will  be  essentially  elastic  as  long 

as  the  shear  modulus  is  maintained  as  the  unloading  modulus  within 
the  cycling  load  range. 

Two  other  triaxial  experiments  were  conducted  where  the 
sampies  were  consolidated  under  a  confining  stress  of  400  psi, 
loaded  vertically  in  the  drained  state  to  630  psi  and  then  further 
loaded  cycled  between  575  psi  and  690  psi  in  both  the  drained  and 
undrained  states.  A  comparison  of  the  results  is  shown  in  Fig.  23, 
in  which  anticipated  responses  were  determined. 
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Fig.  16  Triaxial  Response,  Ranch  Sand  Model,  No  Pore  Fluid 
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Fig.  21  Triaxial  Response,  Ranch  Sand  Model, 
Consolidated  Undrained  Test 
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Fig.  22  Triaxial  Response,  Ranch  Sand  Model,  Effective  Stress  Invariants 
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4.0  SUMMARY 

Numerical  results  presented  herein  fall  into  two  catego¬ 
ries.  The  first  set  of  data  presented  are  concerned  with  determin¬ 
ing  the  adequacy  and  characteristics  of  the  numerical  solutions  for 
combined  stress-seepage.  These  results  indicate  that  the  formula¬ 
tion  and  the  associated  computer  code  developed  to  treat  these  prob¬ 
lems  are  complete  and  debugged.  The  final  set  of  data  is  concerned 
with  attempting  to  evaluate  the  adequacy  of  some  nonlinear  soil  con¬ 
stitutive  models  in  predicting  soil  response  to  load.  The  simple 
Coulomb-Mohr  model  is  clearly  inadequate  except  for  some  simple 
problems  where  strength  alone  is  of  interest,  However,  for  those 
problems  where  stress-history  is  significant,  the  Coulomb-Mohr  model 
must  be  judged  inadequate  except  possibly  to  judge  gross:  strength 
behavior. 

The  McCormick-Ranch  model  (or  types  similar  to  this)  are 
of  course  a  significant  improvement  since  they  will  at  least  repro¬ 
duce  some  known  experimental  responses.  It  can  be  anticipated  that 
they  would  be  adequate  for  various  static  problems  of  interest  or 
for  those  with  only  one  or  two  load  unload  cycles.  However,  these 
models  must  be  judged  inadequate  in  predicting  responses  under 
cyclic  loadings  such  as  those  encountered  in  seismic  problems. 
Further  experimental  data  must  be  developed  for  loading  situations 
with  many  cycles. 
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appendix  a 


formulation  of  SYSTEM  EQUATIONS 
In  the  following  presentation,  the  analysis  will  be  car- 

”16d  £°rth  f°r  3  typiCal  lament  °£  the  free-f ield  mesh.  The  dis¬ 
placement  field  for  the  element  is  assumed  to  be  linear  and  the 
displacement  of  any  point  within  the  element  can  be  -.rittet  as 


*  {%)'{(>} 


(A.  1) 


wnere  (U,w)  are  the  horizontal  and  vertica]  Hie  i 

ana  Vertlcal  displacement  component 

and  W  ,  are  each  a  set  of  arbitrary  coefficients,  with  the 
number  of  coefficients  equal  to  the  number  of  element  vertices  to 
provide  the  proper  number  of  degrees  of  freedom  for  the  element. 

The  vector  {^|  is  formed  by  a  proper  set  of  element  func¬ 
tions  and  depend  upon  the  element  type  being  considered.  For  a 
typical  triangular  element  (Fig.  A.l)  this  vector  is 


while  for  a  typical  rectangular  element  (Fig.  a. 2) 


(A. 2) 


(A. 3) 
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For  a  general  quadrilateral  element  (Fig.  A. 3)  the  vector  V  spec¬ 
ified  by  equation  A. 3  is  used  in  the  transformed  coordinate  system. 
By  substituting  the  coordinates  of  the  nodes  into  equation  A.l, 
the  coefficients  ^  ^  can  be  replaced  as  unknowns  by  the  node 

point  displacement  components,  or 

(A. 4) 

urCr,})  ’  M 

This  simple  displacement  function  assumed  for  the  element  allows 
for  determining  any  interior  displacement  in  terms  of  the  nodal 
displacements  and  ensures  that  the  displacements  between  any  two 
adjacent  elements  will  be  continuous  for  any  arbitrary  specifica¬ 
tion  of  nodal  displacements.  Higher  order  element  formulations 
are  also  available  to  satisfy  the  above  criteria. 

The  strains  developed  at  any  point  within  the  element 
can  be  determined  from  the  strain  displacement  relations  for  the 
particular  configuration,  or 

+ 

where  {£T|  is  the  strain  vector  with  components 

The  superscript  T  in  equation  A. 6  indicates  total  strains. 


(A. 5) 


(A. 6) 
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For  the  combined  stress/pore  pressure  problem,  the  intergranular 
stress  are  related  to  the  pore  pressure  by 


fa]  =  {*]  ' 


(A. 7) 


where  fa|  are  the  total  stress  in  the  body  defined  by 


l<rM<rri  <r.  ,  c-j ,  r,j} 


(A. 8) 


^5"}  are  the  effective  or  intergranular  stresses  and  IT  is  the  pore 
pressure.  The  vector  is  defined  as  ^1,1,1,  o"j  .  The  effec¬ 
tive  stresses  are  related  to  the  strains  through  the  general  stress 
strain  relations 


(A. 9) 


where  are  defined  as  the  nonlinear  components  of  the  total 

strain  and  [c^  is  the  usual  elastic  stress  strain  matrix  which,  for 
example,  can  be  defined  for  the  axisymmetric  problem  by 


Lc]  *  I 
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and  £  1-2*0 
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where  E  is  Young's  modulus  and  })  is  Poisson's  ratio.  For  relatively 
simple  material  models  (such  as  The  Mises  or  Coulomb-Mohr  plastic 
models)  ,  the  nonlinear  strains  represent  the  nonrecoverable  or 
plastic  strain  components.  For  more  complicated  material  models, 
the  vector  represents  a  fictitious  set  of  strains  required  to 

yield  the  proper  stresses. 

To  satisfy  equilibrium  conditions  at  the  element  nodes 
with  the  total  stress  field  within  the  element,  the  usual  virtual 
work  principal  is  used.  The  internal  work  performed  by  the  stresses 
on  a  virtual  total  strain  field  is  defined  by 

Wi  *  ^  (A.H) 

where  the  integral  is  taken  over  the  element  volume.  The  corre¬ 
sponding  external  work  performed  by  forces  applied  at  the  nodes  is 

(S'We  1  (A. 12) 

where  are  the  horizontal  force  components  at  each  node  and 

[r^]  are  the  corresponding  vertical  force  components.  Equating 
the  internal  and  external  work  expressions  and  making  us  of  the 
definitions  previously  described,  the  force  components  that  must 
be  applied  at  each  node  point  to  maintain  equilibrium  with  the 
total  stress  within  an  element  are 


I  i 


!  L  +  Lfeu^H  -  {r“.}  - 
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(A. 13) 


matrices  L uur  1 '  etc.  are  the  usual  elastic  stiffness 

matrices  and  are  defined  by 


Uu]'[  LBcVtcU^lW 


(A. 14) 


where  the  subscripts  U,j)  take  on  the  values  of  (u,w) .  The  terms 

re  iresent  the  correction  forces  to  account  for  mate- 
rial  nonlinearity  and  are  defined  by 


’  I i\i 


(A. 15) 


Where  again  the  subscript  (i,  takes  on  the  vaiues  (u,w).  The  terms 

IX] '  (T<Y1  represent  the  effects  of  pore  pressure  on  the  equilib- 
rium  equations  and  are  defined  by 


Wl  •  I  TT  Ml}  IV 


(A. 16) 


In  the  computer  program  developed,  the  pore  pressure  variation  is 
assumed  to  be  a  linear  one  over  the  element  or 


(A. 17) 


where  the  vector  represents  the  nodal  point  pore  pressures. 
Substituting  equation  A. 17  into  A. 16  then  yields 


where 


ten  'tww 
an  •  £  an'wuyai  * 


(A. 18) 


(A. 19) 


where  the  subscript  l  represents  both  the  u  and  w  directions. 

To  relate  the  pore  presusres  to  the  node  point  displace¬ 
ments  (or  velocities) ,  seepage  effects  are  considered.  The  seepage 
equations  are  obtained  by  minimizing  the  functional  (Ref.  8). 


(A. 20) 


The  permeability  components  |*rr,  are  related  to  the  prin¬ 

cipal  permeability  coefficients  by 


krr-  k,  We  4 
kr4  -s  (te*,-  k,') 

*  kvcosxd  +  hx 


(A. 21) 


where  (k^fe^)  are  the  principal  coefficients  in  two  orthogonal 
directions  and  ©  is  the  angle  from  the  \  -direction  to  the  f -direction 
positive  m  the  clockwise  sense. 

For  a  particular  element,  the  variation  of  pore  pressure 


over  the  element  is  assumed  to  be  linear  or 


i 
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Trtr.jv  'A,22) 

Substituting  equation  A. 22  into  A. 20,  the  functional  can  be  written 


A*iW'ltOls3LiS3H-{'«},LfTi-} 


where  the  matrix  is  defined  by 


(A. 23) 


(A. 24) 


The  vector  ^ j]  is  defined  as  the  volume  integrals 

H  ’  Iv  ^  & 


(A. 25) 


where  Q  is  the  volume  decrease  per  unit  volume  per  unit  time. 

The  solution  to  the  seepage  problem  is  obtained  by  mini 
mizing  the  functional  with  respect  to  the  nodal  pressures  leading 
to  the  conservation  equation 


M'WMH  '  wtj) 


or  to  further  condense  notation 


(A. 26) 


j 


(A. 27) 


The  vector  j  can  be  evaluated  by  considering  the  volume 
compression  of  the  solids  plus  that  of  the  water,  or 

Q  *  '  ^(ST}  4-  g  T  (A. 28) 

where  E  v is  the  effective  bulk  modulus  of  water  and  th'  vector 
^4T}is  the  total  strain  rate  vector.  The  strain  rates  are  related 
to  the  nodal  velocities  by 

[e7]  *  W  r  (A.  29) 

Substituting  the  above  into  the  conservation  equation  leads  to 


IhIV*}  *'  -  fij'i  ir(  +  L-:*l 


(A. 30) 


where 


